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Abstract The hard X-ray modulation telescope mission HXMT is mainly devoted to per- 
forming an all-sky survey at 1 kcV - 250 keV with both high sensitivity and high spatial 
resolution. The observed data reduction as well as the image reconstruction for HXMT 
can be achieved by direct demodulation method (DDM) iteratively. However the original 
DDM is computationally too expensive for multi-dimensional data with high resolution 
to employ for HXMT data. In this article we propose an accelerated direct demodulation 
method adapted for data from HXMT. Simulations are also presented. 

Key words: methods: data analysis — methods: numerical — techniques: image pro- 
cessing — instrumentation: detectors — instrumentation: high angular resolution 

1 INTRODUCTION 
1.1 Hard X-ray modulation telescope 

The hard X-ray modulation telescope (HXMT) is a l ow orbit space telescope runs at a circular orbit 
with a 43° inclinatioin at the altitude of 550 km E l2007l:lLuetalil2010l) . There are three detectors 
on board, including High energy X-ray detector (HE), medium energy X-ray detector (ME), and low 
energy X-ray detector (LE). HE is the most important one, which consists of 18 modules and each 
module is composed o f an HE collimator, a phoswich scintillation detector and readout electronics ( 10 

l2007l:lHanetalll2010h . 

The HE collimator is used to define the field of view (FOV) of the detector dHan et al.ll2010t) . There 
are two types of HE collimators differ in their FOVs: 15 of them are with 5.7° x 1° FOVs, and 3 of them 
are with 5.7° x 5.7° FOVs. Their optical axes are parallel but the directions of their cross-sections are 
different. 

There are two imaging observation modes designed for the science goal, the all-sky survey mode as 
well as the deep imaging observations of selected sky regions. The all- sky su rvey is performed through 
the orbiting of the satellite as well as the precession of i ts orbit p l ane( 0, 120071) . The orbit period is about 



95.5 min while the precession rate is —5.45° per dav dLu et all 1201^ 1 The deep imaging observations 
will be performed by pointed observations with pointing directions distributed uniformly in the region 
and the pointed observations will be performed by progressive scanning. 
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1.2 Modulation 

According to lLi&Wul (Il994 . an observation process can be modeled as 

J p{uj, x)f{x)dx = d{uj), (1) 

where f{x) is the intensity distribution of the object (i.e., the image), is the observed data mod- 
ulated by integral kernel (modulation function) p{uj,x). In practice the observed data is recorded as 
discrete trunks (3-D), maps (2-D) or series (1-D). We have 

N 

dk=Y,Pk,^f^,'^k = l,■■■ ,M, (2) 

1=1 

or 

d = Pf, (3) 

where d is an Af x 1 column vector (data vector) in data space {d} representing all possible observed 
data records, / is an iV x 1 column vector (object vector) in image space {/} representing all possible 
images, thus P is an AI x N matrix, i.e., the kernel matrix in {d} x {/}. 

1.3 Challenges in image reconstruction for HXMT 

Generally speaking, the direct demodulation method (DDM) is designed to solve the demodulation 
problem, i.e., fi nding a best obje ct f{x) to satisfy Eg. [T] whil e the kernel p (a;, x) and the ob served data 
d{ijj) are known jLi & Wuiri993l) . Richardson-Lucy iteration dRichardsoni 1 1 972t ILucvI 1 1 9741) is used for 
DDM, hence the numerical evaluation of modulations is the most time-consuming part. The computa- 
tional complexity increases rapidly as the size of the data increases in multi-dimensional modulation 
evaluation, as a result it is not feasible to reconstruct the all-sky image from HXMT observed data by 
the original DD M. Hence acceleration s of the method are required. 

According to lShen & Zhoul (|2007*) modulations can be accelerated through fast fourier transforms. 
It's implied that the modulation is shift-invariant in this way the modulation equation is reduced into a 
convolution equation. However because of the topology of spherical surface the point spread function 
(PSF) must satisfy certain requirement so that a convolution can be defined there. 

The FOVs of HXMT detectors are not circularly symmetric, i.e., contributions to a detector from 
the sources in its FOV depend upon not only the radial distances from the sources to the center of the 
FOV, but also the directions. In addition, during the all-sky survey the FOV of each detector travels on 
a spherical surface instead of a plane. As a result, the FOVs of each detector around different positions 
on the celestial sphere are never parallel with each other For example, the path of an FOV during the 
first phase of the all-sky survey is shown in Fig [T] Therefore neither fourier transform nor spherical 
harmonics can be used to reduce the modulation of kernel function of HXMT observation and images 
on the celestial sphere. 

In the following sections of this article we present a pixelization and tessellation scheme of data 
defined on spherical surface as well as a method to accelerate the numerical evaluation of modulation 
on such data thus the DDM for HXMT. 

2 IMAGE RECONSTRUCTION FOR HXMT 

2.1 Pixelization and tessellation of data on spherical surface 

Discretization is always the first step in numerical analysis. In image reconstruction pixelization 
is a more specific term for this step. Pixelization of data defined on a plane dr aws little atten- 
tion. But pixelization of data on spherical surface has attracted much m ore interest jTegmarkl 1 19961 : 
ICrittenden & Turokl [19981: iDoroshkevich et all l2005t iGorski et al.ll2005h . Researches on this problem 
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Scanning path in all-sky survey, Phase I 
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Fig. 1 Scanning path of the FOV of a detector of HXMT during the first phase of its all-sky 
survey. We exaggerate the precession rate of the satellite by 20 times to seperate adjacent 
scanning circles in this diagram. 



are driven by applications, e.g., CMB data analysis, since due to the topological nature of a spherical 
surface, an theoretically i deal pixelization scheme that would work for all cases does not exist for data 
defined on such a surface (iTegmarki [19961 iGorski et al.Ll2005l) . 

Image reconstructions for both all-sky survey and observations for selected sky regions can be 
conquered region by region locally since given a single point of the observed data, only objects within 
the current FOV are relevant while given a single point on the unknown image, only local observed data 
is relevant. Therefore all-sky pixelization is not necessary. We need only a local pixelization scheme 
adequate for a sky region larger than the FOV of the telescope, which is 5.7° x 5.7°. 

Virtues of a pixelization scheme in this work may include; 

- All pixels are of equal size, which speeds up and simphfies the evaluation of numerical integration, 
the elementary operation of modulations. 

- Geodesic between two points p ^ {i^, iy) (ix and iy are indices) and q — i'^,iy as well as azimuth 
of one point with respect to the other are both shift-invariant. That is, 

d[iix,iy), (i'x, iy)] = d[{ix + jx, iy + jy), (i'x + jx, iy + jy)], yix,iy, i'xJyJxJy, (4) 

where d{p, q) is the geodesic distance between p and q, i.e., the length of the minor arc on great 
circle from p to q, and, 

Ct[{ix,iy), (i'xJy)] = (5) 

where a{p, q) is the angle where the geodesic arc from p to g crosses the meridian containing 
p. Therefore all pixels should be of the same shape. If the PSF of the telescope is circularly 
symmetric, i.e., the PSF can be expressed as a function of distance from the center of the FOV, the 
modulation then degenerates to convolution. 
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- 2-D orthogonal pixel indexing. For example the pixels are indexed with two indices and iy in 
Eq.|4]and Eq.|5] where it suggests that the pixel indexing consists with a 2-D Cartesian coordinate 
system, which is the default indexing way in numerical analysis of 2-D image. Therefore the indices 
ix and iy serve as dimensionless orthogonal coordinate variables in the discrete model where the 
image value: 

where Sx and Sy are the sampling intervals along x-axis and y-axis respectively, while U^.i^ is 
the sampled image value on {ix,iy) and I{x, y) is the image value defined on {x, y) in the orginal 
continuous model. 

Equidistant cylindrical projection (ECP) method is commonly used in geophysics and climate mod- 
eling. Equidistant pixels on the surface of a cylinder c are projected to the surface of a inscribed sphere 
s of the cylinder c towards their symmetry axis. Adjacent pixels on the spherical surface are either of 
the same right ascension (R.A.) or of the same declination (Dec). Pixels on the same parallel (or merid- 
ian) are spaced on adjacent meridians (or parallels) uniformly. Pixels with higher latitudes (closer to 
poles) have smaller sizes and greater distorsions in contrast with pixels with lower latitudes (closer to 
the equator). 

HEALPix bv lGorski et al.l (l2005h has recently become a standard structure for spherical data anal- 
ysis especially for CMB experiments. Although the equal areas of HEALPix pixels is vital for spherical 
harmonics transforms, however, the shapes of pixels are different thus the shift-invariance is violated. 
Therefore HEALPix is not adapted for the accelerated image reconstruction method we proposed for 
HXMT data in this article. 

2.7.7 Quadrilateral projection based pixelization 

A pixelization scheme is designed for HXMT data. We use the scheme introduced here to build a grid 
of pixels on a specific sky region. There are several concepts we need to specify: 

1 . A pixel on a sphere is an elementary area on the sphere inside its boundary and around its center. 

2. The boundary of a pixel on a sphere is defined by 4 vetices of the pixel and 4 geodesies of the 
sphere between each adjacent pairs of them. 

3. The position of the center of a pixel on a sphere is fixed to indicate the area inside the boundary of 
the pixel. For example, given a extreme case where the 4 vetices lie on the same great circle of the 
sphere thus the great circle itself yields the boundary, the two part of the spherical surface divided 
by the great circle are identical. Therefore the common sense that a pixel is a small portion of the 
whole image does not help to define the pixel. So a point inside the pixel must be given to indicate 
the area inside. 

We start with generating pixels with equal areas and the same shape in a square on a tangent plane 
p of the unit sphere. This area is divided into Np ^ N x N square pixels. Then the centers of all the 
pixels are projected from the plane p to the surface of the unit sphere s. The projection can be either 
radial or parallel. In the radial mode, each pixel center is projected towards the center of the sphere s, 
while in the parallel mode it is projected perpendicularly to the plane p. 

Once all the centers of the pixels are fixed on the sphere, find the geodesies between each center 
and its 4 nearest neighbours. Next draw a great-circle arc across the middle point of each geodesic 
perpendicularly. The 4 points where each perpendicular great-circle arc intersects another two arcs are 
the vertices and the 4 perpendicular great-circle arcs make the boundary. With the center indicating the 
inside area of the boundary, a pixel is defined on the sphere around each center 

Pixels on the spherical surface projected far from the square center have greater distorsions in con- 
trast with those projected from the central area of the square. Obviously the pixel distorsion is indepen- 
dent of the the position of the pixel on the sphere but only depends on its original position in the plane, 
more specifically, the distance from the square center on the plane, therefore we can always use square 
small enough to make the distorsion negligible. 
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We use the term tessella to refer to a set of pixels on the spherical surface, which are projected 
from all pixels within a square on a plane. Hence it suggests that we divide the problem pixelization 
of a whole spherical surface into two problems, pixelization of any small region of the spherical 
surface and tessellation of the whole spherical surface using tessellae of pixels. Generally speaking 
there should not be overlaps or gaps between adjacent tessellae, however, since we will not perform 
numerical evaluations on different tessellae in the same time, overlaps only make evaluations on pixels 
of the overlaps redundant. 

Since all tessellae have the same area and shape, we can generate a set of pixels (i.e. a tessella) 
around the null position of the spherical surface and rotate this initial tessella to a series of positions to 
cover the whole sphere. In this way we can pixelize any small region of a spherical surface by the initial 
tessella and its rotation expressed by a quaternion (Appendix [All. 

As shown in Fig. |2] we define that the null position of the surface of a unit sphere s is (1, 0, 0) and 
the plane pis x — 1. The center of this square is projected to (1, 0, 0) perpendicularly to p. 




Fig. 2 Plane projection around null position of the surface of a unit sphere. The unit sphere 
s is centered on (0, 0, 0), the origin of the xyz 3-D Cartesian coordinate system. The null 
position is on (1, 0, 0). The plane p is parallel to the plane yOz. There are 4 pixels on the 
plane. They are either radially or parallelly projected to s. 



As shown in Fig. [3] Fig. |4] and Fig.|5] the pixel distorsions are negligible in a small region close to 
the equator, for all the three projection-based pixelization schemes. The distorsions become significant 
while the region expends. Although in low latitude region the distorsion of ECP scheme is suppressed 
better than those based on plane projection, soon we find that in high latitude regions the ECP scheme 
suffers from severe distorsion, even in a small region. But obviously the distorsions of plane projection 
schemes are independent of latitudes. 

See statistics on shapes and areas of pixels on the sphere in Fig. |6] 
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Fig. 3 Pixel distorsion of equidistance cylindrical projection. Left: 11.25° x 11.25°, centering 
at 0° latitude. Middle: 30° x 30°, centering at 0° latitude. Right: 11.25° x 11.25°, centering 
at 60° latitude. 



Fig. 4 Pixel distorsion of radial plane projection. Left: 11.25° 
latitude. Middle: 30° x 30°, centering at 0° latitude. Right: 11. 25=^ 
latitude. 



X 11.25°, centering at 0° 
x 11.25°, centering at 60° 



Fig. 5 Pixel distorsion of parallel plane projection. Left: 11.25° x 11.25°, centering at 0° 
latitude. Middle: 30° x 30°, centering at 0° latitude. Right: 11.25° x 11.25°, centering at 60° 
latitude. 



2.7.2 Tessellation scheme 

The size of the tesseUa is cr x cr, where 

0-= M = 1,2,.... 

2M 

For the radial projection the size of the square on the tangent plane is 2tan^^ | x 2tan~^ ^. 
tessellation of the unit sphere is implemented in 3 steps. 



Accelerated DDM for HXMT 



7 



X 10" 




relative distance angle, degrees 



Fig. 6 Pixel statistics. 1 024 x 1 024 pixels are projected to 11.25° x 11.25° spherical surface 
radially. Top: relative pixel area, proportion of the average of areas of all the pixels; bottom 
left: relative pixel distance, proportion of the average of distances between adjacent pixels; 
bottom right: angle between adjacent sides of the boundary of each pixel. 



The first step is tessellation of the prime meridian of the sphere. Rotate the initial tessella from the 
null position towards the north pole as well as the south pole at intervals of cr, as shown in Fig.|2] Then 
the prime meridian is covered by 2M + 1 tessellae. We call each tessella of those the meridianal tessella 
(Fig.Q. 

The second step is tessellation of the equatorial area of the sphere. Rotate the initial tessella around 
z-axis also at intervals of a so that the equator of the sphere is covered by 4M equatorial tessellae, as 
shown in Fig. [8] 

The last step is tessellation of the remaining area of the sphere. Rotate each meridianal tessella 
except the initial one and the two covering the polar caps around z-axis at intervals of 
where 0^ is the latitude of the center of each tessella, i.e.. 



6^ = ± , m = 1, .... M - 1. 

See the tessellation along a given latitude in Fig.|9] 

Finally the number of tessellae covering the whole sphere is 



(8) 



A/-1 r 



Nf = AM + 2- 



tt/ arctan( 



tan 



4M 



COS + Sin tan 3^ 



(9) 
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Fig. 7 Meridianal tessellae on the northern celestial hemisphere 




Fig. 8 Tessellation of equatorial area and polar area 



where [a:] is the ceiling function denoting the smallest integer no less than x, Vx G 5R. See the tessella- 
tion of the celestial sphere in Fig.fTOland Fig.fTTI 

The quaternion representing the rotation upon each tessella is 

= cos cos — + sm — sm — i - sm — cos — j + cos — sm — k, (10) 
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Fig. 9 Tessellae along a given latitude 




Fig. 10 Tessellation of the northern celestial hemisphere 



where and 9^ represent the position of the center of the tessella on the sphere, in longitude and 
latitude respectively. The corresponding rotation matrix is 

(cos9m cos(j)w — sin(/)^ — sin6'u, cos^^.X 
cos sin (/)^ COS0U, — sin 6'^, sin 0.^, . (11) 
sin^^ cos6'u, / 

Given a pixel rij on the initial tessella, where i and j are the indices of each pixel, the corresponding 
pixel on the tessella around {(pw, dw) is then R^fy^^g^ri^j. 

2.1.3 Pixelization of observed data 

The original observed data from each detector module is recorded as a TOA (time of arrival) sequence 
of x-ray photons. It is required to convert the original observed data from time series into pixel-wise 
format, i.e., pixelization of the observed data, for following processes. Alongside with the scientific data 
there is also engineering data being recorded such as the attitude of the spacecraft etc. By interpolations 
the attitudes of the spacecraft can be determined at any given point during the observation time, i.e., for 
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Fig. 11 Tessellation of the celestial sphere 



each detected photon the attitude of the spacecraft can be assigned. From the attitude of the spacecraft 
we can calculate the position on the celestial sphere where each collimator is pointing. For example, 
in the equatorial coordinate system, given a specific detector module, each detected photon from this 
module comes with not only its TOA but also the R.A. and Dec. of a position on the celestial sphere 
the collimator is pointing at, i.e., the projection on the celestial sphere along the optical axis of the 
collimator. The projection on the celestial sphere is then called the position of the detected photon. 

For each photon first we find the tessallae that cover its position on the celestial sphere. Then for 
each tessella that covers its position, we find the pixel that takes up this position. Because the boundary 
of each tessella as well as for each pixel is defined, the problem of finding whether a specific point on the 
celestial sphere lies inside or outside a tessella or a pixel is classified as a point-in-square-on-sphere 
problem. See in Appendix IB] our approach to solve this problem. 

Once the position of each detected photon is mapped to pixels and tessellae we defined on the 
celestial sphere, the total number of detected photons can be calculated for each pixel. Taking the 
time interval between the time of arrival of a photon and the next one on the same detector as the 
exposure time of the former photon, the exposure time of all the photons mapped to each pixel together 
yields the total exposure time on this pixel. Divide the total number of photons by the total exposure 
time on each pixel, and we have the counting rate for each pixel, i.e., the observed data in pixel-wise 
format. 
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2.2 Acceleration of modulation 

2.2. 1 Status parameters of the coUimated detector 



Evaluating numerical modulation is ubiquitous in both simulation of observed data and image recon- 
struction with direct demodulation method. In HXMT observation the observed data is a function of the 
given collimator status, which includes the collimator identifier c and its orientation. 

To describe the orientation of a given collimator, we use a unit vector along its optical axis as well 
as another unit vector along a fixed axis perpendicular to its optical axis (e.g., long edge of one of its 
slice). We call the unit vector along its optical axis and the later one pointing vector and position vector 
(imitating the position angle parameter of a telescope) respectively. Then we define a null status of a 
given collimator where its pointing vector points at the null position of the current coordinate system 
(along .T-axis) and its position vector points at (0, 0, 1) (along z-axis). We use three status parameters 4), 
9 and ^ to represent the current status of the given colhmator. First rotate the given collimator around 
2:-axis from its null position by </> then incline it from the xOy plane (e.g. the equator) by 6, and finally 
rotate it around its own optical axis by V' i-c, the position angle. In this way the colUmator can be rotated 
to any possible orientation. 

We use a quaternion q to formulate the rotation parameterized by (j), 9 and if) on the pointing vector 
and position vector of a given colUmator as 

6 6 

9i = cos I + sin |k, (12) 

9 9 9 

= cos - + sine/) sin -i — coS(/)sin -j, (13) 

^3 = cos — + cos (j) cos ^ sin — i + sin cp cos 6 sin — j + sin 6 sin — k, (14) 
q = 939291. (15) 

where q^, and q^ are auxiUary quaternions, and q rotates the colUmator from its null status ( to 
^ I ^ . The corresponding rotation matrix is 

(cos 9 cos (f) — sin (j) cos ip — sin 9 cos (j) sin ip sin sin ip — sin 9 cos (j) cos ip \ 
cos^sint/) cos (/) cos i/l" — sin sin sin — cos </) sin f/; — sin sin </> cos -0 . (16) 
sin^ cos ^ sin ^ cos ^ cos V' / 

Given the status parameters (p, 6 and tp, the corresponding pointing vector and position vector are 

/l\ / cos 9 cos (f)\ 
P4>,0,^ = R<P,e,->P = <^os9tim(f) , (17) 



, / \ sin ( 



and 



sin (j) sin ip — sin 9 cos 4> cos -0 

0'4,,e,ii> = R<t>,e,ip I I = I — cos ^ sin — sin sin ^ cos | , (18) 

cos 6 cos V' 

respectively. 

We call the vector ^ e ^ the status vector of a given collimator, where <p, 9 are the two pointing pa- 
rameters and ^p is the position angle parameter. The status parameters are coordinate-system-dependent. 
Null status vector ^ o ^ in different coordinate systems stands for different orientations of the same 
coUimator. 
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For example, if the equatorial coordinate system is adopted, ip ^ means its position vector and 
pointing vector are on the same longitude of the celestial sphere, i.e., the R.A. of both its position 
vector and pointing vector are equal to each other The angle between the two vectors are always f by 
definition. So in this case the Dec. of its position vector differs by 90° from the Dec. of its pointing 
vector. 



2.2.2 PSF and numerical evaluation of the modulation kernel function 

In the equatorial coordinate system an image is a function of positions on the celestial sphere denoted 
by R.A. and Dec. 

Given the observed data d{(t>, 6, ijj, c) and the image /((/>, 6), the modulation in Eq.[T]yields 

9, V, c) = / p{^, 9, i;, c, (t>\ 0')f[<P\ 0')An, (19) 



Jn 

where 9, ijj, c, (p' , 9') is the modulation kernel function represents the response of collimator c with 

/ \ / COS 9 cos (j)' \ 

status I e ) to a unit point object at position I COS 9' sin (j)' ) , and respresents the solid angle. 

^ sine' ' 

The PSF of collimator c is defined through its modulation kernel function as 

P(0,0,c) 0,0, c, 0,0), (20) 

i.e., the response of collimator c to a unit point source located at the null position of the celestial sphere 
when the collimator is pointing at {(\),9) and its position angle parameter remains zero. In practice the 
PSF is measured during the calibration of the detector(,Han et al., , 2010) , estimated from observations, 
or predicted theoretically. In any case the PSF of a given imaging detector is determined by its intrinsic 
properties and characterizes its response to a unit object. The response to the unit object is determined 
only by the position of the unit object relative to the collimator rather than their absolute positions in any 
coordinate system. So the modulation kernel function in a specific coordinate system can be evaluated 
through the PSF P((/), 9, c) by coordinate transforms. 

To evaluate the modulation kernel function 9, ip, c, 0', 9') we find a coordinate system where 

the coordinate of the unit object is o ^ > i-e., the null position of the new coordinate system and the 

position angle parameter of the collimator yields zero. Once we find the pointing parameters $ and 
of the collimator in the new coordinate system, the modulation kernel function can be evaluated 
immediately. 

This is implemented by the following steps: 

1 . Rotate the current coordinate system 5*0 to the first auxiliary coordinate system Si where the unit 
object lies on o ^ For example, we can first rotate the orginal coordinate Sq around its z-axis by 

(/>' then arount its y-axis by —9' . 

2. Find the status parameters (f>i, 9i, and iIji in Si. The coordinate of a vector in the rotated coordinate 
system 5*1 is equivalent to rotating the vector inversely, i.e., rotating it around the z-axis of 5*0 by 
—cf)' then rotating around y-axis of 5*0 by 9' . We formulate the rotation by quaternions as: 

91 = cos y - sin y k, (21) 

9' 9' 

92 = cos y + sm yj, (22) 

q = 929i' (23) 
where qi and ^2 ^re auxiliary quaternions. The corresponding rotation matrix is 

cos 9' cos (f)' cos 9' sin 0' sin 9' ^ 
= ( -sin0' COS0' |. (24) 

■ sin 9' sin 6' cos 9' , 
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Given that the pointing vector and position vector of the collimator are g ^ and a^.e.ip in So as 
stated in Eq.[T7]and Eq.[T8]respecively, the two vectors of the collimator are 

Pi = Re'.,4>'P^.0,^, (25) 

and 

o-i = Re' ,4>'Ci'it>,e,ii)- (26) 

3. Rotate Si around its x-axis by —a so that the position angle parameter of the given collimator in 
the rotated coordinate system 52 is zero. Given the pointing vector and the position vector of the 
collimator in ^2 are ™d 02 respectively, we have (02 x P2) • k = 0, i.e., the cross-product of 
the position vector and the pointing vector lies on the xOy plane of 5*2. Therefore we have 

a = Arg{Rg,^^,{a^^e,i, x P^^e,-4,) ' K Re',4,'{a^,e,i, x P^^e,4,) ' j)' (27) 
where Arg(a;, y) is the principal argument of a complex number with x + yi. 

4. Find the pointing parameters (t>2 and 02 of the collimator in S2 and use them to evaluate the PSF, 
which is equivalent to the modulation kernel function in 6*2. Given the rotation matrix 

/I \ 
= COS a — sin a , (28) 
\Q sin a cos a / 

where a is solved in Eq.|271 the pointing vector P2 of the collimator is 

P2 ^ Ra^R<P',e'Pci,fi,r (29) 
Hence we can find the pointing parameters (/)2 and 62 from 

02 =Arg(p2 • j,P2 • i)' (30) 
62 =arcsin(p2 • k). (31) 
Finally the modulation kernel function is evaluated as 

p{(^,e,^,c,c^',e')^p{cj)2.e2,c). (32) 

2.2.3 Numerical evaluation of modulation 

Since the response of a detector of HXMT to an object depends upon the position of the object relative 
to the detector only, we can always rotate both the detector and the object to the initial tessella of the 
celestial sphere to calculate the modulation kernel function, if the tessella is larger than the FOV of the 
detector, i.e., 

P(P<?i,e,V" a^.e.V" = PiR^l,e^P4>,e,i,^ c, -R^^,e„P0',e')' (33) 

where p^ g ^ and a,t>fi,^ represent the pointing vector and the position vector of the collimator which 
are related to its status (0, 6*, V'), P^'.e' represents the point {(/)', 9'), so p{p^ g ,^,a^^g^^,c,p^, g,) is 
equivalent to 9, ip, c, 0', 9'). Rotation matrix rotates the initial tessella to {(f>w, (^w)- 

The first-order approximation of the modulation kernel function around (f) ~ Q, 9 — Q, (f)' — and 
9' = (i.e., on the initial tessella) is 

p[c^, 9, V', c, </.', 9') = {R^P){(j) ~c^',9~ 9', c) + iV((/.2) + N[9^) + N{(j)''^) + N[9'% (34) 

where we put a rotation matrix to the left of the PSF P to define a new PSF {R^P) which is rotated 
from the original one. 

Taking the pixelization of observed data into account, the corresponding discrete modulation equa- 
tion yields 

AT 

^i;j,c — ^^{R^ijP^i — i',j—j\cfi\j'-! (35) 

where ipij is the position angle parameter of the collimator c on {i, j) and Pi.j^c is the value of the PSF 
of the collimator on 
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2.2.4 Position angle cluster analysis and approximation 

Because of the topology of the spherical surface, the position angle parameter of each collimator varies 
during the all-sky survey and the variance depends upon the position on the celestial sphere. Given a 
specific scanning scheme for the all-sky survey, e.g., the satellite orbits along the path shown in Fig.[T] 
and its roll angle is fixed on a constant value such as 0°, —30° or 30°, the position angle parameter of a 
collimator is either t/^c ~ 43° or ipc + 43° (V'c is the angle between the position vector of the collimator 
and the instantaneous velocity of the satellite, i.e., the tangent vector of the scanning trajectory) approx- 
imately in the initial tessella, but varies in a wide range significantly in some other tesselae, as shown in 
Fig. El 
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Fig. 12 Position angle variance and distribution. Top left: position angle variance along scan- 
ning path of all-sky survey in equatorial area; top right: position angle distribution in equa- 
torial area; bottom left: positon angle variance along scanning path of all-sky survey in high 
latitude area; bottom right: position angle distribution in high latitude area. 



On the pixel grid of a small sky region we can rotate the PSF around itself and shift it along two 
orthogonal directions to evaluate the modulation kernel function approximately, as shown in Equation 
[35] Because of the position angles variance as well as the lack of circular symmetry of the PSF the 
modulation kernel function is shift-variant. We can approximate the modulation kernel function with 
the sum of a finite series of functions so that the position angle parameter of the collimator status is 
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fixed for each function, as 



/ \ ' ^^^^ 

fc=i ^ / " 

where ^^^^K is rounded to the nearest integer to evaluate the Kronecker delta 5^, equivalent to 

the Dirac delta function S{il; — 

We use cluster analysis to assign the position angle parameters of all the observed data in a given 
sky region into groups so that the difference between each position angle parameter of a group and 
the clustering center of the group is less than a pre-determined upper limit e according to the required 
precision. Therefore the modulation kernel function is shift-invariant approximately if all of its possible 
position angle parameters are in the same group. So we rewrite Equation|35]by expanding the right-hand 
side: 

N 
i' , ?' 

(37) 

= ^^(-RVfe-P)i-«'j-j',c/*'j'(5(|V'fc - Aj\ < e) 

fc=l i' J' 

where 

Finally we approximate the modulation by the sum of a finite series of convolution as 

K 

dij,c = ^[(-Rv'.Pc) * fkjSil^k - V'^.jI < e), (39) 

where [{R^^Pc) * f] represents the convolution between the rotated PSF of collimator c and the image 
/. Then we can employ FFT algorithms to accelerate the modulation. 

3 SIMULATED DATA AND RECONSTRUCTED IMAGES 
3.1 PSF of the HXMT main detector 

The PSF of collimators as well as the HE detector are shown in Fig.fTSl 



3.2 Simulated all-sky survey data and reconstruction 

The 7-year INTEGRAL all-sky survey catalos dRrivonos et al. I I2OO7I) is used as input to simulate HXMT 
data, as shown in Fig. [14] 

Two small regions are selected for simulated observations as well as reconstructions. One of them 
is in equatorial area while the other is in high latitude area. The azimuth variance in high latitude area is 
more significant thus the simulated observation as well as the reconstruction consume more computation 
time. 

The model images, observed data as well as the reconstructed images are shown in Fig. [15] 

To test the resolving ability of DD method for HXMT all-sky survey data we set up a point source 

with 1 000 mCrab flux as well as a uniform background so that the significance of the source is ka, 

where a is the standard deviation of the fluctuation of the background. 

We simulated the observed data and the reconstructions for all the set-ups, as shown in Fig. [16] 
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Fig. 13 PSFs of HXMT main detector. Left: PSF of collimator with 5.7° x 1° FOV; middle: 
PSF of collimator with 5.7° x 5.7° FOV; right: PSF of the high energe detector, overlaid. 



INTEGRAL 7-year survey 




R.A., rad 



Fig. 14 Sources in 7-year INTEGRAL all-sky survey catalog. Upper box: a tessella in the 
equatorial area. The R.A. and Dec. of its center are —80° and 0° respectively. Lower box: a 
tessella in the high latitude area. The R.A. and Dec. of its center are —101° and —37°. 



4 CONCLUSIONS 

In this article we present the approximation of the kernel function for image reconstruction from HXMT 
data and the pixelization scheme optimized for this accelerated DD method. The pixels are nearly 
squares with the same size to make up a grid suitable for convolution. However due to the spherical 
topology it is impossible to tessellate the spherical surface without gap or overlap with perfect spherical 
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squares of the same size. The error of the sizes of the pixels are propagated to the results of numerical 
integrations. The distortions of the pixels cause the error of distance between any two pixels, which 
are propagated to the numerical convolutions eventually. These errors will not be accumulated or mag- 
nified t hrough iteration s because each iteration of DD method is based on the modulation equation 
directlv (lLi&WuLll993h . 

With simulated data we demonstrate the method we proposed works for both equatorial data and 
data observed in high latitude area. We can obtain a reconstructed image from the given data in several 
minutes on ordinary desktop PCs. 

Since we are focused on the acceleration of the DD method, some other topics are not discussed, 
such as noise suppression, treatment of faint and extended sources, which are also necessary for image 
reconstruction for HXMT data. For example, in Fig. [15] we adjust the SNR of the observed images 
so that images that with the model images can be achieved although a number of the sources in the 
INTEGRAL 7-year catalog are actually too faint for a single phase of the HXMT all-sky survey. In 
addition we use known backgrounds as background constraints in the DD method, which is usually too 
idealized to realize in practice. 

The acceleration technique we proposed in this article is necessary and sufficient for our future 
work. It serves as a tunable underlying engine of the flexible DD method. 

Appendix A: QUATERNION AND RELATED OPERATIONS 

In 3-D Cartesian coordinate system a quaternion is defined as 

q = a + bi + ci + dk, (A.l) 

where a is its scalar part and 6i + cj + dk is its vector part. 

A quaternion q = a + 6i + cj + dk is a unit quaternion if and only if + 6^ + + = 1. 
A unit quaternion is used to formulate a rotation performed on a rigid body in 3-D space. Given a unit 
quaternion 

q = a + M + cj + rfk = cos — + sin — (cos 9 cos (f>i + cos 9 sin ip] + sin 9k), (A. 2) 

/ cos 9 cos (j) \ 

q then indicates the right-handed rotation around vector ( cos e sin <p = cos 9 cos 0i+cos 9 sin (p] +sin 9k 

V sin e / 

by angle a. 

The conjugate of quaternion q = a + &i + cj + dk is defined as 

q = a — &i — cj — dk. (A. 3) 

If the quaternion q is unit then its conjugate q is equivalent to its inverse q~^. 

Given the equation of basis elements i, j and k of quaternions as well as 3-D vectors 

i^ = = k^ = ijk = -1, (A.4) 

the multiplication of quaternions and 3-D vectors can be calculated distributively. 

Let q and t> be a unit quaternion and an arbitrary 3-D vector respectively, hence the multiplication 
qvq^ yields the vector we obtain on rotating v as q indicates. 

Appendix B: POINT-IN-SQUARE-ON-SPHERE PROBLEM AND 
SPHERICAL-RAY-CASTING ALGORITHM 

The point-in-square-on-sphere (PISOS) problem asks whether a given point on the spherical surface 
lies inside, outside or on the boundary of a spherical square. A spherical square is a regular spherical 
quadrilateral which has four equal sides and four equal angles. Each side is the arc of great circle passes 
two adjacent vertices while each angle is formed by the tangents of two adjacent sides. 
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The PISOS problem arises when we try to decide which observed data or pixel lies in which sky 
region. This problem is a special case of PIP (Point-in-polygon) problem in computational geometry, 
which can be tackled by the ray-casting algorithm. According to the algorithm, one can find whether a 
given point is inside or outside a polygon by testing how many times a ray, starting from any known 
point inside the polygon and going towards the given point, intersects the edges of the polygon. If the 
number of the intersections is odd the point is outside while it is inside if even. 

We designed a more specific method, the spherical-ray-casting algorithm, to solve the PISOS prob- 
lem. Assuming the spherical square is smaller than a semispherical surface, the PISOS problem reduces 
into the following problems: 

1. Find intersections of two given great circles on a unit sphere. Either great circle is defined by two 
points on the unit spherical surface. 

This problem is reduced into finding the normal vector of a plane defined by two points on a unit 
spherical surface as well as the orgin. For example, let p^, and p^ be points on the unit 

sphere -\- + = 1, where m is the normal vector of the plane p^Opj 'ind n2 of p^Op^, 
i.e., rii = I p^^p^ i and n2 = jf^^fj- If 9 is the intersection of plane PiOp2 and p^Op^, vector q 
is perpendicular to both rii and n2 thus any vector on the plane n\On2, therefore q is the normal 
vector of this plane, i.e., q = 

r ' ' Im Xn2 1 

Let 6 1, (f)i and 6*2, ^2 be the altitudes and azimuthal angles of two points on the unit sphere while 
On and (j)n be those of their normal vector. We have 

cos 6n COS (pn COS 6*1 COS (pi + COS On Sin (/)„ COS Oi sin 01 + sin On sin 6*1 = 
cos On cos (pn COS O2 COS 02 + COS On sin (/)„ COS O2 sin 02 + sin On sin 02 = 

hence. 



D = cos^ ^1 sin^ O2 + sin^ Oi cos^ 62 — "2 cos cos 62 sin sin 02 cos(0i — 02) 
cos 61 COS 62 sin(02 — 0i) 

, (B.2) 



tan On = 



sm 0n = 



COS0„ = 



D 

sin Oi cos O2 cos 02 — cos di sin 62 cos 0i 
D 

cos Oi sin O2 sin 01 — sin Oi cos O2 sin 02 



D 



where D is an auxiliary quantity. 
2. Find whether a given minor arc of great circle contains a point on the same great circle. Let pg and 
Pi be the endpoints of an minor arc of great circle and let g be a point on the great circle. The point 
q is on the given arc if and only if the sum of geodesic distances from q to Pq and to Pi equals to 
the length of the arc. 
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Fig. 15 Top left: Model image in equatorial area. Top right: Model image in high latitude 
area. Middle left: Observed image in equatorial area, with poisson noise. Middle right: 
Observed image in high latitude area, with poisson noise. Bottom left: Reconstructed im- 
age using accelerated DD method, with 100 000 iterations (about 3 minutes). Bottom right: 
Reconstructed image using accelerated DD method, with 100 000 iterations (about 10 min- 
utes). 
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Fig. 16 The upper 4 images: Observed data of point sources with 20(T, 5(t, 3<t and la- 
significances. The SNRs are 11.1 dB, 5.1 dB, 2.88 dB and -2.07 dB respectively. The 
lower 4 contour diagrams: Reconstructed images of the point sources. 



